System and method for interpolating seismic data by matching pursuit in fourier transform

ABSTRACT

Methods and systems for interpolating seismic data estimate a first frequency spectrum associated with the seismic data using a matching point (MP) algorithm having a pre-computed term. A second frequency spectrum is also estimate for the seismic data, using the MP algorithm, an anti-aliasing mask and the estimated first regular frequency spectrum. From the second frequency spectrum, the interpolated seismic data can be determined by performing an inverse fast Fourier Transform thereon.

PRIORITY INFORMATION

The present application claims priority under 35 U.S.C. §119(e) to U.S. Provisional Patent Application Ser. No. 61/805,259, filed 26 Mar. 2013, the entire contents of which are expressly incorporated herein by reference.

TECHNICAL FIELD

The present embodiments relate generally to land and marine seismic exploration systems and methods and, more specifically, to systems and methods for interpolating seismic data using matching pursuit algorithms.

BACKGROUND

Seismic waves generated artificially have been used for more than 50 years to perform imaging of geological layers. During seismic exploration operations, vibrator equipment (also known as a “source”) generates a seismic signal that propagates in the form of a wave that is reflected at interfaces of geological layers. These reflected waves are received by geophones, hydrophones or, more generally, receivers, which convert the displacement of the ground resulting from the propagation of the waves into an electrical signal which is recorded. Analysis of the arrival times and amplitudes of these waves make it possible to construct a representation of the geological layers on which the waves are reflected.

FIG. 1 depicts schematically marine seismic exploration system (marine system) 100 for transmitting and receiving seismic waves intended for seismic exploration in a marine environment. Marine system 100 comprises a source 118 on a streamer or cable 116 a, pulled from ship or boat 102, on the surface 106 of ocean 108 (or other water mass, such as a large lake or river). Source 118 is operable to generate a seismic signal. Marine system 100 further includes a set of receivers 120 (e.g., hydrophones) for receiving a seismic signal and converting it into an electrical signal, also located on streamer 116 b, and marine seismic data recording/processing system 126 for recording and processing the electrical signals generated by receivers 120. Streamers 116 can also include birds 122 for guiding and maintaining position of streamers 116. Source 118, receivers 120 can be intermixed on one or more streamers 116, in any order. FIG. 1 depicts source 118 as a single source but it should be understood that the source may be composed of several sources, as is well known to persons skilled in the art. Also part of system 100 are antennas 124 that can be used to transmit information and controls between ships 102, land bases, birds 122 (some birds 122 can also include antennas 122) and other devices.

In operation, source 118 is operated so as to generate a seismic signal. This signal propagates through water 108, in the form of transmitted waves 124 that generate reflected waves 126 when they reach an interface 110 between two layers 108 (ocean) and 112 (a geological layer, in this case, the ocean floor; it can also be appreciated by those of skill in the art that sometimes the transmitted waves 124 travel upwards instead of downwards, and when they reach the interface between atmosphere/air 104 and ocean 108 (i.e., at ocean surface 108) downward reflected waves 126 can also be generated (not shown); these are known by those of skill in the art as “ghosts”). Each receiver 120 receives one or more reflected waves 126 and converts them into an electrical signal. Marine system 100 intends to image the subsurface regions 112 to determine the presence, or not, of hydrocarbon deposit 114.

FIG. 2 depicts a land seismic exploration system 200 for transmitting and receiving seismic waves intended for seismic exploration in a land environment. Land system 200 comprises a source 202 consisting of a vibrator operable to generate a seismic signal, a set of receivers 204 (e.g., geophones) for receiving a seismic signal and converting it into an electrical signal and land seismic data recording/processing system 226 for recording and processing the electrical signals generated by receivers 204. Land system 200 can further include antennas 124 for communications between vehicles 226, receivers 204, and land seismic data recording/processing system 226.

Source 202, receivers 204 and land seismic data recording/processing system 224 (located on vehicle 226) are positioned on the surface of ground 208. FIG. 2 depicts source 202 as a single vibrator but it should be understood that the source may be composed of several vibrators, as is well known to persons skilled in the art. In operation, source 202 is operated so as to generate a seismic signal. This signal propagates firstly on the surface of the ground, in the form of surface waves 210, and secondly in the subsoil, in the form of transmitted waves 212 that generate reflected waves 214 when they reach an interface 220 between two geological layers. Each receiver 204 receives both surface wave 210 and reflected wave(s) 214 and converts them into an electrical signal, which signal thus includes a component associated with reflected wave 214 and another component associated with surface wave 210. Since land system 200 intends to image the subsurface regions 216 and 218, including hydrocarbon deposit 214, the component in the electrical signal associated with surface wave 210 is undesirable and should be filtered out.

It is known to those of skill in the art that seismic processing techniques often assume data are regularly sampled in space, but acquisition methods, particularly in marine environments, rarely achieve this in practice. Seismic data which exhibits, for example, large gaps or has irregular sampling, needs interpolation in order to fully benefit from further processing.

Seismic data interpolation methods can be roughly divided into those which are model-driven from the intrinsic wavefield physics, and those which are data-driven and exploit patterns in the data. In the former category, wave-equations use wavefield continuation methods to integrate new traces at nominated positions, but often require knowledge of the wavefield propagation velocity. In the latter category, signal processing based approaches assume that seismic traces contain locally coherent events, whose patterns allow energy to be projected from neighboring input traces to construct new traces. Examples of signal processing approaches include spatial Prediction Error Filter (PEF), which are designed from input traces and able to predict missing data traces (see, Spitz, S., 1991, “Seismic Trace Interpolation in the f-x Domain,” Geophysics, 56, 785-794), and a whole class of algorithms that construct regularized seismic data from input data transformed into some intermediate domain. The transform is usually Fourier, but it can also be a Radon or a curvelet transform.

Recently, Fourier transform methods for seismic data interpolation have become popular because they are relatively easily extended to higher dimensions. As seismic data is generally well sampled in the time direction, such data can be transformed stably into the f-x domain. For each temporal frequency f; a 3D interpolation algorithm attempts to estimate the spatial Fourier spectrum (i.e. the f-k_(x)-k_(y) domain). The increased data from additional dimensions (e.g. offset and azimuth) stabilizes the transform estimation (see, Trad, D., 2009, “Five-Dimensional Interpolation Recovering from Acquisition Constraints,” Geophysics, 74, V123-V132). However, this technique is subject to several constraints. For example, its inverse Fourier transform should exactly or approximately fit with available f-x data. Additionally, as described in “Seismic Trace Interpolation in the Fourier Transform Domain,” Gulunay, N., Geophysics, 68, 355-369 (2003), the cyclic property of the discrete Fourier transform imposes a predefined shape on the spectrum.

Another example of Fourier transform based interpolation is the minimum weighted norm interpolation (MWNI) method (see, Liu, B., et al., 2004, “Minimum Weighted Norm Interpolation of Seismic Records,” Geophysics, 69, 1560-1568), which minimizes an objective function that includes two terms: i) the difference between interpolated and observed data at input coordinates, and ii) a weighted spectral norm which constrains the shape of the reconstructed spectrum. Still another Fourier transform interpolator is the anti-leakage Fourier transform (ALFT) (see, Xu, S., et al., 2005, “Anti-leakage Fourier Transform for Seismic Data Regularization,” Geophysics, 70, V87-V95). The Xu method is considered by those of skill in the art to be simple and intuitive: given input samples in the f-x domain, estimate their discrete Fourier transform in the f-k domain. The largest Fourier coefficient is selected and subtracted from the input data. In subsequent iterations, successive maximum components are subtracted until the norm of the residual is negligible. The motivation for the subtraction is to attenuate “frequency leakage” of the most energetic Fourier coefficient to other coefficients potentially attributable to irregularly sampled data. This iterative process is able to recover a sparse spectrum that, when evaluated at sampling points, approximates regularly sampled data. Several of these methods are discussed in greater detail below. However; there is a fundamental problem with each of the aforementioned processes: they rely on the common assumption that sparsely sampled data can be adequately represented by a any Fourier components.

This common assumption, however; is not universally true in the world of seismic data acquisition. As can be appreciated by those of skill in the art, seismic processing techniques often assume data are regularly sampled in space, but acquisition methods, particularly in marine environments, rarely achieve this in practice. Furthermore, aliasing problems can occur when reconstructing seismic data. If a seismic event is too deep, for example, in relation to other events (“event” or “events” meaning in this context reflection sources), then the high frequency values of the trace data become abased, i.e., folded over or mixed up with lower frequency data.

Accordingly, for seismic data that exhibits large gaps and/or irregular sampling, it would be desirable to provide methods, modes and systems for interpolation that overcomes the problems in the conventional techniques described above.

SUMMARY

An aspect of the embodiments is to substantially solve at least one or more of the problems and/or disadvantages discussed above, and to provide at least one or more of the advantages described below. It is therefore a general aspect of the embodiments to provide an interpolation technique for seismic that will obviate or minimize problems of the type previously described.

According to a first aspect of the embodiments, a method for generating interpolated seismic data includes the steps of estimating a first frequency spectrum for seismic data using a matching point (MP) algorithm having a pre-computed term, estimating a second frequency spectrum of the seismic data, using the MP algorithm, an anti-aliasing mask and the estimated first frequency spectrum; and determining the interpolated seismic data by performing an inverse fast Fourier Transform on the second estimated frequency spectrum.

According to a second aspect of the embodiments, a system for generating interpolated seismic data includes at least one processor configured to estimate a first frequency spectrum for seismic data using a matching point (MP) algorithm having a pre-computed term, to estimate a second frequency spectrum of the seismic data, using the MP algorithm, an anti-aliasing mask and the estimated first frequency spectrum; and to determining the interpolated seismic data by performing an inverse fast Fourier Transform on the second estimated frequency spectrum.

According to a third aspect of the embodiments, a method for generating 5D interpolated seismic data using a fast matching pursuit (FMP) algorithm includes the steps of: obtaining a block of Np irregular trace data, and defining N_(k) to be a desired number of regular trace data; estimating a conjugate of a Gram Matrix G* with size N_(p)×N_(k), using the N_(p) block of irregular traces; transforming the obtained irregular trace data from a t-x domain to an f-x domain; estimating a first set of frequency slices f₁ in an f-k domain of the transformed trace data using a fast MP algorithm that uses G*, wherein f₁ is for low frequency values only; estimating an anti-aliasing mask M based on the first set of frequency slices f, estimating a second set of frequency slices f₂ using the fast MP algorithm using G* and the anti-aliasing mask M, wherein f₂ includes all frequencies up to and including a Nyquist frequency of the transformed trace data; and determining N_(k) interpolated trace data by reverse transforming f₂ from the f-k domain to the t-x domain using an inverse 5D fast Fourier Transform.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and other aspects of the embodiments will become apparent and more readily appreciated from the following description of the embodiments with reference to the following figures, wherein like reference numerals refer to like parts throughout the various figures unless otherwise specified, and wherein:

FIG. 1 illustrates a side view of a marine seismic exploration system for use in an underwater seismic gathering process;

FIG. 2 illustrates a side view of a land seismic exploration system for use in a land seismic gathering process;

FIG. 3 illustrates a general method for seismic exploration according to an embodiment;

FIG. 4 illustrates the decay of the norm of the residual R^(m)ƒ, of an anti-leakage Fourier transform algorithm and a fast matching pursuit algorithm according to an embodiment, plotted in logarithmic scale;

FIG. 5 illustrates a comparison of the Fourier spectrums produced by an anti-leakage Fourier transform algorithm and a fast matching pursuit algorithm according to an embodiment, plotted in logarithmic scale;

FIGS. 6A-6C illustrate examples of synthetic trace data interpolated using an embodiment of the fast matching pursuit algorithm;

FIGS. 7A and 7B illustrate examples of actual trace data interpolated using an embodiment of the fast matching pursuit algorithm;

FIG. 8 illustrates a flow chart of a method for 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment;

FIG. 9 illustrates a method for determining a set of regular frequency slices using a Gram Matrix, G*, according to an embodiment;

FIG. 10 illustrates a method for determining an anti-aliasing mask according to an embodiment;

FIG. 11A illustrates a flow chart of a method for 5D anti-aliasing matching pursuit seismic data interpolation according to another embodiment; and

FIG. 11B illustrates a flow chart of a generalized method for interpolating seismic data according to another embodiment;

FIGS. 12A-12D illustrate examples of synthetic trace data interpolated using 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment;

FIGS. 13A and 13B illustrate examples of actual trace data interpolated using 5D anti-aliasing matching pursuit seismic data interpolation algorithm according to an embodiment;

FIGS. 14A and 14B illustrate stacks prior to, and after, use of the 5D anti-aliasing matching pursuit seismic data interpolation algorithm, respectively, according to an embodiment; and

FIG. 15 illustrates a marine seismic data acquisition system suitable for use to implement a method for 5D anti-aliasing matching pursuit seismic data interpolation algorithm according to an embodiment;

DETAILED DESCRIPTION

The embodiments are described more fully hereinafter with reference to the accompanying drawings, in which embodiments of the novel concept are shown. In the drawings, the size and relative sizes of layers and regions may be exaggerated for clarity. Like numbers refer to like elements throughout. The embodiments may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be complete, and will convey the scope of the concepts therein to those skilled in the art. It will be apparent to one skilled in the art, however, that at least some embodiments may be practiced without one or more of specific details described herein. In other instances, well-known components or methods are not described in details or are presented in simple block diagram format in order to avoid unnecessarily obscuring the embodiments. The scope of the embodiments is therefore defined by the appended claims. The following embodiments are discussed, for simplicity, with regard to the terminology and structure of a marine seismic exploration system. However, the embodiments to be discussed next are not limited to these systems but may be applied to other, e.g., land seismic exploration systems that are affected by problems of seismic data interpolation.

Reference throughout the specification to “one embodiment” or “an embodiment” means that a particular feature, structure, or characteristic described in connection with an embodiment is included in at least one embodiment of the embodiments. Thus, the appearance of the phrases “in one embodiment” on “in an embodiment” in various places throughout the specification is not necessarily referring to the same embodiment. Further, the particular feature, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

According to the embodiments to be described below, matching pursuit seismic data interpolation algorithm methods overcome problems associated with the conventional techniques described in the Background section. For example, some of these interpolation embodiments make use of the Gram Matrix to reduce the computational complexity significantly over conventional interpolation methods. More specifically, and according to an embodiment, the number of computations can be reduced from Np×Nk operations to only Nk operations, wherein Np is the number of input traces, and Nk is the desired number of output traces. Additionally, some of the described embodiments include an anti-aliasing feature that further enhances the quality of the interpolated results and reduces (or eliminates) reliance upon the assumption that sparsely sampled data can be adequately represented by a few Fourier components.

Prior to discussing these interpolation embodiments in more detail, it may be useful to consider the overall seismic exploration process in general for context. As generally discussed above, one purpose of seismic exploration is to render the most accurate graphical representation possible of specific portions of the Earth's subsurface geologic structure, e.g., using the seismic data which is collected as described above with respect to FIGS. 1 and 2. The images produced allow exploration companies to accurately and cost-effectively evaluate a promising target (prospect) for its oil and gas yielding potential (e.g., hydrocarbon deposits 114). FIG. 3 illustrates a generalized method for seismic exploration that includes both the acquisition of the seismic data described above, and the subsequent processing of that seismic data to form such images. In FIG. 3, the overall process 300 is broken down into five steps, although one could of course characterize seismic exploration in a number of different ways. Step 302 references the initial positioning of the survey equipment in the geographic area of interest (GAI) and the preparation to begin surveying the GAI in a manner which is precise and repeatable. Seismic waves are generated by the afore-described sources or vibrators (step 304), and data recording is performed on the reflected, scattered and surface waves by the receivers (step 306).

In step 308, processing of the raw, recorded seismic data occurs. Data processing generally involves numerous processes intended to, for example, remove noise and unwanted reflections from the recorded data, and involves a significant amount of computer processing resources, including the storage of vast amounts of data, and multiple processors or computers running in parallel. In particular, for the embodiments discussed below, such processing includes interpolation to “fill-in” missing and/or incomplete and/or noisy recorded data. Such data processing can be performed on site, back at a data processing center, or some combination thereof. Finally, in step 310, data interpretation occurs and the results can be displayed or generated as printed images, sometimes in two-dimensional form, more often now in three dimensional form. Four dimensional data presentations (i.e., a sequence of 3D plots or graphs over time) are also possible outputs, when needed to track the effects of, for example, extraction of hydrocarbons from a previously surveyed deposit.

With this context in mind, and as briefly discussed above, there are numerous methods available for performing seismic data interpolation to fill in gaps within an acquired seismic data set. One such method can be generally referred to as matching pursuit. Described below is a new, fast matching pursuit (FMP) algorithm according to an embodiment, followed by a comparison between the FMP algorithm with a conventional, anti-leakage Fourier Transform (ALFP) method which is also used to perform seismic data interpolation.

To illustrate a method of regularization of non-uniform sampled functions by matching pursuit, interpolation problems associated with a single frequency slice of 2D seismic data are first considered. Throughout this discussion, bold letters imply column vectors and bold capital letters represent two-dimensional arrays. A continuous complex-valued function ƒ(x) can be defined in the range from 0 to 1 and sampled at a set of N_(p) points, as follows:

xl, 0≦xl<1, 1≦l≦N _(p).

Any continuous function ƒ(x) can be decomposed into an infinite sum of complex exponentials that are periodic in the range [0, 1]. It can be assumed that the function is band-limited and can be expressed as a sum of N_(k) frequency components from −N_(k)/2 to +N_(k)/2−1. These complex exponentials can be evaluated exactly by the integral Fourier series formula only if ƒ(x) is known for all x. The values {circumflex over (ƒ)}(k), −N_(k)/2≦k<N_(k)/2 are called. Fourier coefficients of ƒ(x). If it is desired to find {circumflex over (ƒ)}(k) so that the sum of the N_(k) complex exponential functions at x_(l) are equal to ƒ(x_(l)), this can be accomplished by computing:

$\begin{matrix} {{f\left( x_{l} \right)} = {\sum\limits_{k = {{- N_{k}}/2}}^{\frac{N_{k}}{2} - 1}\; {{\hat{f}(k)}{^{2\pi \; \; {kx}_{l}}.}}}} & (1) \end{matrix}$

Equation (1) expresses the N_(p) values of ƒ(x_(l)) as the sum of N_(k) Fourier functions e^(2πkx) ^(l) . In vector notation, samples ƒ(x_(l)) can be written as a column vector ƒ of length N_(p); the N_(k) frequency components are written as column vector {circumflex over (ƒ)}. The exponential wave function φ_(k) evaluated at x_(l) is denoted by k, with elements:

${\frac{1}{\sqrt{N_{p}}}^{2\pi \; \; {kx}_{l}}},$

wherein the term:

$\frac{1}{\sqrt{N_{p}}}$

normalizes the vector φ_(k) to have a norm of 1. The matrix φ of size N_(p)×N_(k) is created by N_(k) column vectors φ_(k). Equation (1) can be expressed in matrix form as:—

ƒ=√{square root over (N _(p))}φ{circumflex over (ƒ)}  (2).

Equation (2) is an undetermined system of equations when N_(k)>N_(p). The use of the Fourier method in seismic interpolation assumes that the spectrum is sparse, which means there are relatively few dominant wavenumbers, most of the spectral values being near zero.

One way to find a sparse solution to an undetermined system of equations is by using the matching pursuit (MP) algorithm. Those of skill in the art will appreciate that MP is a so called ‘greedy algorithm’ that optimizes the approximation of ƒ by selecting one vector φ_(k) at each iterative step. Effectively, the MP algorithm searches for the index p₀ of the set of φ_(k) that best correlates with ƒ, as:

|(ƒ,φ_(p) ₀ )|≧|(ƒ,φ_(k))|,−N _(k)/2≦k<N _(k)/2  (3).

The number of wave vectors N_(k) is determined by the required output grid for interpolated data output. The vector ƒ can be written as the sum of its orthogonal projections onto the vector φ_(k) and the residual R¹ƒ:

R ¹ƒ=ƒ−(ƒ,φ_(p) ₀ )φ_(p) ₀   (4).

The MP algorithm continues by selecting another vector φ_(p) ₁ that best estimates the residual vector R¹ƒ. MP attempts to find a sparse approximate solution {circumflex over (ƒ)} to Equation (2) using the method discussed below. The sparse approximate solution {circumflex over (ƒ)} is obtained when the residual vector R¹ƒ is estimated to be less than the threshold ε, as discussed below. The threshold value ε is normally calculated from the seismic data and a user parameter. According to one embodiment, at the start of the process, the program reads the user parameter, for example 0.001, and then multiplies it with the average of the absolute value of a random set of 100 input seismic traces. This value is used as epsilon. If the user chooses a low threshold (e.g., less than 0.001), then the program will take a longer time to run, and vice versa. The MP method for finding a sparse approximate solution {circumflex over (ƒ)} to Equation (2) begins with choosing a threshold ε, setting m=0, {circumflex over (ƒ)}=0, R^(m)ƒ=ƒ (step 1). The MP algorithm then finds (step 2) the index p_(m) of the vector φ_(p) _(m) that best correlates with the residual R^(m)ƒ by estimating:

{tilde over (ƒ)}^(m) =Φ*R ^(m)ƒ  (5),

where {tilde over (ƒ)}^(m) is a vector of size N_(k), and its p_(m) element is:

{tilde over (ƒ)}(p _(m))=

R ^(m)ƒ,φ_(p) _(m)

  (5A).

In a third step the MP algorithm then increments {tilde over (ƒ)}^(m)(p_(m)) by

$\frac{1}{\sqrt{N_{p}}}{{{\overset{\sim}{f}}^{m}\left( p_{m} \right)}.}$

The MP algorithm then subtracts (in step 4) from the existing residual value R^(m) the new incremented value of {tilde over (ƒ)}^(m)(p_(m)):

R ^(m+1) ƒ=R ^(m)ƒ−{tilde over (ƒ)}^(m)(p _(m))φ_(p) _(m)   (5B),

and determines the norm of the new residual R^(m+1)ƒ. In the final step of the MP algorithm, the MP algorithm determines whether the norm of the residual R^(m+1)ƒ is less than ε, the predetermined threshold, and if so, the MP algorithm has determined the index p₀ of the set of φ_(k) that best correlates with ƒ, as described above. If, however, the norm of the residual R^(m+1)ƒ is equal to or greater than the predetermined threshold ε, the MP algorithm increments m by one, and returns to the step of obtaining the index p_(m) of the vector φ_(p) _(m) .

Those of skill in the art will further appreciate that the most expensive phase in the MP algorithm, in terms of processing bandwidth, is the step where matrix multiplication is required. Such matrix multiplication results in a complexity of each iteration on the order of N_(p)×N_(k). According to an embodiment, however, if the equation in step 4 is multiplied by Φ* and using of Equation (5), then it can be determined that:

{tilde over (ƒ)}^(m+1)={tilde over (ƒ)}^(m)−{tilde over (ƒ)}^(m)(p _(m))Φ*φ_(p) _(m)   (6).

If the vector Φ*φ_(p) _(m) is pre-computed, then {tilde over (ƒ)}^(m+1) can be estimated from {tilde over (ƒ)}^(m) in only N_(k) operations. As those of skill in the art can appreciate this modification to the matching pursuit algorithm can save a significant amount of processing operations per iteration. As such, the matching pursuit algorithm which is modified as described above with respect to Equation (6) can now be referred to as the fast matching pursuit (FMP) algorithm and represents one aspect of an interpolation technique for seismic data according to these embodiments.

Pre-computation of the vector Φ*φ_(p) _(m) can be performed, according to an embodiment, by pre-computing the so-called Gram matrix. By definition, the Gram matrix G of the set of N_(k) vectors φ_(k) is an N_(k) by N_(k) matrix whose elements g_(i,j)=(φ_(i),φ_(j)). It can be observed by those of skill in the art that the vector Φ*φ_(p) _(m) is the p_(m) column of the conjugate of the Gram matrix, G*, which can be determined as shown in Equation 6(A):

G*=Φ*Φ  (6A).

In the context of seismic interpolation processing, the cost of estimating and storing G* is essentially negligible because G* is a Hermitian Toeplitz matrix with diagonal elements equal to 1. Therefore, only the N_(k) elements of the first row or column need to be computed. More importantly, as Fourier seismic interpolation is run in the f-x domain, the set of basis function φ_(k) stays the same for each ƒ, thus the pre-computed matrix G* can be used for all frequencies f.

In order to better understand the improvement provided by FMP techniques according to these embodiments, relative to conventional interpolation techniques, a comparison will be provided with respect to the conventional anti-leakage Fourier transform (ALFT). Numerical testing of the FMP and the ALFT algorithms was performed using a function having two different frequencies:

ƒ(x _(l))=5 sin(25.6x _(l))+2.5 sin(128x _(l)+1)  (7).

As those of skill in the art can appreciate, the is the same test function used by Xu et al., as discussed in his article published in 2005, “Anti-Leakage Fourier Transform for Seismic Data Regularization,” Geophysics, 70, V87-V95, which the interested reader may also consult to obtain details regarding the ALFT, generally. The function in equation (7) was sampled randomly at N_(p)=64 points in the range [0,1] and the two algorithms (FMP and ALTF) were used to recover N_(k)=128 frequencies.

Using the results of this numerical testing, FIG. 4 illustrates the decay of the norm of the residual R^(m)ƒ, plotted in log scale, for each of the FMP and ALFT algorithms. Therein, Line A illustrates the decay in residual value for the ALFT algorithm, and line B illustrates the decay in residual value for the FMP algorithm according to an embodiment. In order to reduce the effect of randomness of the sampling points, the plotted decay rates in FIG. 4 represent averaged results after 1000 runs. As can be seen in FIG. 4, the FMP algorithm has a much faster decay rate compared to the ALFT. In general, the ALFT needs to run twice the number of iterations compared to the FMP algorithm with the same stopping criterion. According to an embodiment, the FMP algorithm looks for the best basis function to project the residual onto. Therefore, the norm of the residual function is reduced by the maximum possible amount at each iteration. For the ALFT algorithm, there is no guarantee that the norm of the residual will decrease, and experiments have shown that the norm of the residual can actually increase over a number of iterations when the input data contains noise. Since almost all real world data has noise, this clearly demonstrates an inherent deficiency in the ALFT algorithm.

It is known that the ALFT algorithm has a good approximation of the spectrum at each iteration step, and therefore those of skill in the art could assert that it might be able to produce a better overall spectrum than the FMP algorithm. However, numerical experiments using both the FMP and ALFT algorithms clearly show that the spectra produced by the two algorithms are nearly identical, as FIG. 5 illustrates. Therefore, those of skill in the art can appreciate that the FMP algorithm is advantageous in that it requires less processing than the ALFT to produce similar spectra.

To illustrate the application of FMP embodiments to the interpolation of seismic data, FIGS. 6A-6C depict examples of synthetic seismic trace data interpolated using the FMP algorithm according to an embodiment. The synthetic data consists of a set of 32 traces containing three linear events. Sixteen traces are randomly removed and the MP algorithm is used to reconstruct the whole data set. FIG. 6A is the original set of 32 traces; FIG. 6B is the same set of synthetic data with sixteen (50%) randomly removed traces, and FIG. 6C is the result of applying the FMP algorithm according to an embodiment on the synthetic data set of FIG. 6B with half of the trances randomly removed.

A comparison between FIG. 6A (the original data) and FIG. 6C (the interpolated data) shows the results of the MP algorithm as being very good at reconstructing the original data set of traces. Those of skill in the art can appreciate that a small amount of coherent noise is present, which becomes even more evident if more input traces are discarded. This occurs when the spectrum of the original data in the f-k domain cannot be reconstructed perfectly; any error will be expressed as a coherent wave in the t-x domain. This is a specific problem with Fourier methods, demonstrated clearly here on synthetic data. Other problems are the wrap-around effects noticed on the quiet trace on the extreme right hand side. This can be ameliorated by tapering, which is a process known to those of skill in the art.

FIGS. 7A and 7B illustrate an example of actual trace data interpolated using the FMP algorithm according to an embodiment. More specifically, FIGS. 7A and 7B represent the application of a 3D FMP interpolation algorithm according to an embodiment to a real marine dataset in midpoint coordinate space. The nominal inline sampling is 12.5 m, but with only a 25 m crossline. FIG. 7A illustrates one crossline section, and the maximum gaps are 10 to 15 traces wide and are typically present in all crosslines. The interpolation results of FIG. 7B are fairly good, with potential diffraction, un-conforming events, and the amplitude spectrum is well preserved in the interpolated traces. As those of skill in the art can appreciate, it is difficult to visually determine which traces are interpolated and which traces are input data. Those of skill in the art can appreciate that large gaps are convincingly interpolated.

Attention is now directed towards another matching pursuit algorithm for interpolation of seismic data according to an embodiment. FIG. 8 illustrates a flow chart of method 800 for 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment. Prior to discussing method 800 for 5D anti-aliasing matching pursuit seismic data interpolation in detail, a brief discussion of 5D data and interpolation, in general, will be provided.

It will be appreciated by those skilled in the art that seismic data acquisition produces 5 dimensional (5D) data, in which one dimension is time, and the rest of the dimensions are the coordinates of the source and receiver on the surface. Whenever the input data is 5D, such as 3D land or wide-azimuth marine surveys, a 5D interpolation method is preferable as it will use all available data to reconstruct traces in regular source and receiver coordinates. It is further known that many interpolation methods are transform-based, meaning that the interpolated data is generated by inverse transformation. In the attempt to reconstruct missing data, the main expenditure of processing time/effort is in the construction of the transformed representation that best reconstructs the input data based on in inverse transformation.

Attention is now directed towards FIG. 8, and method 800 for 5D anti-aliasing MP seismic data interpolation according to an embodiment, which begins with step 802, wherein seismic waves (marine or land) are transmitted by sources 118, 202. In step 804, reflected waves, some primary and some as a result of multiples, and in the case of land seismic exploration, ground waves, are acquired by receivers 120, 204, processed, and stored in system 126 (marine)/226′ (land), both of which are described in greater detail below.

In method step 806, method 800 obtains a block of N_(p) irregular trace data from the acquired trace data, and defines N_(k) to be a desired number of regular trace data. That is, N_(k) is the desired output of regular trace data of method 800. In method step 808, method 800 uses the irregular N_(p) coordinates to first estimate Φ, the matrix of N_(p)×N_(k), created by N_(k) column vectors φ_(k). the rows indexed by x_(l) having no particular order according to an embodiment, while the columns indexed by four dimensional wavenumber k follow lexicographical order. The column vectors φ_(k) can also be defined as the exponential wave function k, evaluated at x_(l), wherein x_(l) is a sampled point value of the trace data. Then, method 800 can estimate a conjugate of a Gram Matrix G*, which has a size N_(k)×N_(k), using said N_(p) block of irregular traces.

Next, in method step 810, method 800 transforms the N_(p) irregular trace data from a t-x domain to an f-x domain. In step 812, method 800 then estimates a first set of regular frequency slices f₁ using the fast MP algorithm discussed above according to an embodiment, which uses G* and further wherein the first set of regular frequency slices f₁ is for signal frequency values only, i.e., the frequency of the transmitted seismic signals. The final output from step 812 of method 800 puts the data into the f-k domain. According to an embodiment, signal frequencies can be in the range of about 10 Hz to about 70 Hz, though those of skill in the art can appreciate that these frequencies should not be construed in a limiting sense, in that different frequency values can be used and still considered within aspects of the embodiments.

To better understand the complex processing which is being performed in step 812, FIG. 9 illustrates method 900 for determining a set of regular frequency slices using a Gram Matrix, G*according to an embodiment. According to this embodiment, G* is used to calculate the first (regular) frequency slices in the following iterative manner. A “dirty” spectrum is first determined or estimated from the irregular trace data (steps 902-904). In this context, a dirty spectrum is the spectral frequency range of the irregular trace data, and can be denoted as F_(d). According to an embodiment, F_(d) is determined using any of well-known Fourier transform methods to determine the frequency coefficients. The F_(d) coefficients are obtained in the f-x domain, while the original trace data is in the t-x domain. Following transformation using a Fourier transform, method 900 takes the frequency data, d, of each frequency slice f, and multiplies it by Φ, determined in step 808 of method 800, to obtain F_(D)(n) in the f-k domain, according to the following relationship:

F _(D)(n)=Φ×d.  (8)

Then, using an iterative approach, a regular spectrum, F_(r), is determined (steps 908-920). The regular spectrum is the spectrum or spectral frequency range of the desired complete interpolated set of trace data. The regular spectrum value is initially set to 0: F_(r)=0; and a maximum coefficient value of the dirty spectrum is obtained, F_(d-max), this, of course, indicates some spectral frequency of the dirty spectrum. F_(D-max)(n) is placed in F_(R)(n) at the same location as F_(D-max) was located (step 910). The column of G* that corresponds to F_(d-max) is then multiplied by F_(d-max) (which can then be denoted as G*_(max)), and subtracted from F_(d) (step 912). This produces a new set of dirty spectrum (as well as incrementally and iteratively building up the regular spectrum G*), and the process continues by iteratively taking the largest remaining coefficient value F_(d-max) until none are left, producing the complete G*.

According to an embodiment, the use of different columns of G*, as found in step 912 of method 900 includes the use of the FMP algorithm described above. According to a further embodiment, however, the complete MP algorithm can be used instead, albeit many more computations would be required without obtaining significantly better results. As those of skill in the art can appreciate, however, this illustrates that various other means exist for implementing method 900 for determining a set of regular frequency slices using G* and such other means are considered to be included as different aspects of the embodiments described herein.

Returning again to FIG. 8, and method 800, in step 814, method 800 uses the first set of regular frequency slices f₁, determined in method step 812, to estimate an anti-aliasing mask M (AA Mask). According to an embodiment the AA Mask M can be determined as follows. Generally, the AA Mask M is determined using an iterative approach wherein at each iteration the amplitude of the residual spectrum is multiplied with the previously determined (or known) AA Mask M, and the position of the maximum element in the result is determined. The coefficient of the residual spectrum at this position is selected, and the corresponding component of the wavenumber is subtracted from the residual spectrum.

More detail regarding this technique for estimating an AA mask is provided in FIG. 10. Therein, method 1000 begins with step 1002 wherein, for the first iteration, an initial AA Mask, AA Mask₀ is set equal to all 1's, to initialize it. Further in step 1002, method 1000 sets counter m to 1, and determines an iteration counter limit N, where N is set by a user based on prior experience and/or empirical data, as is known to those of skill in the art. Further in step 1002, method 1000 determines threshold, T, according to a random selection of K trace data point. Once the K trace data points are selected, method 1000 transforms them to the f-k domain, and obtains their average, which becomes the threshold, T. K is a number selected by the user, again based on prior experience and/or empirical data. The threshold T is used to determine, in one embodiment, a stopping point for the iterative determination of the AA Mask; another stopping point is the number of iterative determinations performed, i.e., the number of loops. According to another embodiment, the iterative counter limit N determines the number of iterations (or calculations of the AA Mask) that should be performed, to determine the AA Mask.

Further in step 1002, method 1000 sets the factor α=1. According to an embodiment, the factor α can be set by the user to be a number to be between 0 and 1, depending on the amount of noise predicted to be in the trace data. Method 1000 then proceeds, still while in step 1002, to obtain Φ*, the conjugate of Φ, determined in method 800, step 808, and then to use that result to determine the Gram matrix conjugate G* according to Equation (6A), i.e., G*=Φ*Φ.

Method 1000 then proceeds to method step 1004, wherein an initial residual spectrum is determined according to Equation (5), i.e.:

{tilde over (ƒ)}^(m) =Φ*R ^(m)ƒ;

where Φ* is the conjugate of Φ, and the residual spectrum, RS_(M)=R^(m)ƒ. In step 1006, method 1000 then multiplies the previously determined AA Mask_(M−1), in the initial case, AA Mask₀, against the residual spectrum, to determine the residual spectrum mask product, RS-Mask Pro duct, RSMP_(M).

Then, in step 1008, method 1000 finds the maximum value and corresponding position in the RS-Mask Product, Max_Val(RSMP_(M)) and Max_Pos(RSMP_(M)), respectively. In step 1010, method 1000 uses the Max_Pos(RSMP_(M)) as determined in step 1008 to find its corresponding column in G*; method 1000 then obtains the wavenumbers from that column, which are represented as G*_column_(M). In step 1012, a new residual spectrum, RS_(M+1) is calculated according to the following expression:

RS _(M+1) =RS _(M)−[(α)×(Max_(—) Val(RSMP _(M)))×(G*_column_(M))];

where, Max_Val(RSMP_(M)) was determined in step 1008, and G*_column_(M) was determined in step 1010. As discussed above, the factor α, while typically set to 1, could be different, for example, a value between 0 to 1 dependent upon the amount of noise in the trace data, as previously discussed.

Then, method 1000 proceeds to decision step 1014 wherein a determination is made as to whether stopping criteria has been encountered. That is, method 1000 determines whether the subtrahend of [(α)×(Max_Val(RSMP_(M)))×(G*_column_(M))] is greater than the threshold, T. If [(α)×(Max_Val(RSMP_(M)))×(G*_column_(M))] is greater than T (“Yes” path from decision step 1014), method 1000 is complete, and the last determined residual spectrum is the anti-aliasing mask for use in the balance of steps remaining in method 800. If, however, [(α)×(Max_Val(RSMP_(M)))×(G*_column_(M))] is not greater than T (“No” path from decision step 1014), method 1000 first sets AA Mask_(M) equal to RS_(M), and then proceeds to step 1016, wherein the counter m is incremented by 1, and then checked, in step 1018 as to whether iteration counter limit, N, has been met.

That is, if m equals N, then the iterative loop of method 1000 has been determined to have occurred enough times, and the last determined residual spectrum is used as the anti-aliasing mask (“Yes” path from decision step 1018). The iteration counter limit N is determined by a user as discussed above. If, however, m<N (“Yes” path from decision step 1018), then method 1000 proceeds to step 1006, and using the just determined new residual spectrum, and previous AA Mask_(M−1), determines a new AA Mask_(M). The output of method 1000, therefore, is the anti-aliasing mask M used in step 816 of method 800.

While FIG. 10 represents one way to determine an AA Mask according to an embodiment, it will be appreciated that the method of FIG. 8 is not limited thereto, and that other embodiments of FIG. 8 may employ known techniques for determining an anti-aliasing mask.

Returning to FIG. 8, method 800 now uses the whole range of frequencies f, in step 816, to estimate a second set of frequency slices f₂ using the same fast MP algorithm as used in step 812, using G* as determined in method step 808 and the anti-aliasing mask M determined in method step 814, wherein f₂ includes all frequencies up to and including a Nyquist frequency of the transformed trace data. According to an embodiment, the anti-aliasing mask M is used to weight the residual in order to preferentially selected the most energetic k coefficient as observed from method step 810, according to the MP algorithm as discussed by Nguyen, T., et al., in an article entitled “Seismic Interpolation by Optimally Matched Fourier Component,” 2011, SEG International Exposition and 81th Annual Meeting, 3085-3089. According to an embodiment, the most energetic k coefficient is the coefficient with the greatest amplitude following weighting. The second set of regular frequency slices f₂ is the basis for the interpolated data. In this context, the second set of frequency slices f2 refers to all of the frequency slices needed to reconstruct interpolated traces. For example if the traces are sampled at 4 ms, the Nyquist frequency is 125 Hz and the second set of frequency slices f2 include frequency slices from 0 to 125 Hz.

According to an embodiment, the final AA Mask M determined in method step 814 is used to estimate the second set of frequencies, the regular frequency slices, f₂. AA Mask M is used to weight the residual (or dirty spectrum) by multiplication. That is, AA Mask M is multiplied against the value of the residual, which, according to an embodiment, are complex numbers. Then, the coefficient at the position where the weighted value is maximized is selected. As will be appreciated by those skilled in the art, the term ‘frequency slice’ refers to the whole spectrum that the MP algorithm is trying to reconstruct. For example, suppose that a 2D MP algorithm is trying to reconstruct 32 seismic traces of length 100 samples in the t-x domain. The MP algorithm does that by reconstructing an f-k representation of 100 ‘frequency slices’ f, each frequency slice consisting of 32 frequency k. So the MP algorithm will be running 100 times with the antialiasing mask for each frequency slice. Accordingly, to estimate the second set of frequencies, the present embodiment multiplies the antialiasing mask with the current residual frequency, and finds the position that has maximum value for a specific frequency k. This can be accomplished by, for example, calculating, for a selected frequency index n, the equation:

n=argmax_(—) n(Mask*absolute_value(Phi*R̂mf))

In method step 818, method 800 determines the interpolated data by reverse transforming f₂ from the f-k domain to the t-x domain using an inverse 5D Fast Fourier Transform (FFT), the inverse 5D FFT being known to those of skill in the art. In optional step 820, method 800 performs further processing and/or displaying of interpolated data.

According to a further embodiment, method 800 can be implemented as method 1100, shown in FIG. 11A, which also accomplishes data interpolation. Method 1100 for 5D anti-aliasing matching pursuit seismic data interpolation begins with step 1102, wherein seismic waves (marine or land) are transmitted by sources 118 or 202. In step 1104, reflected waves, some primary and some as a result of multiples, and in the case of land seismic exploration, ground waves, are acquired by receivers 120 or 204, processed, and stored in system 100 (marine)/200 (land), both of which are described in greater detail below.

Next, in step 1106, an estimation is made of a regular frequency spectrum for irregular seismic data using the FMP algorithm. Then, in method step 1108, method 1100 determines an anti-aliasing mask M based on the estimated regular frequency spectrum as determined in method step 1106. In method step 1110, method 1100 re-estimates the regular frequency spectrum of the irregular seismic data, using the fast MP algorithm and the anti-aliasing mask M. In method step 1112, method 1100 determines the 5D interpolated data by performing inverse fast Fourier Transform on the re-estimated, regular frequency spectrum. Further processing and/or displaying of interpolated data can then occur in step 1114.

Various alternatives are also contemplated by the embodiments. For example, the estimation of the anti-aliasing mask and the actual MP interpolation described above can be performed independently of one another. For example, the MP interpolation described above can be used with a different anti-aliasing mask, or with a dip mask estimation method. A more generalized embodiment is illustrated in the flow diagram of FIG. 11B. Therein, at step 1120, a first frequency spectrum for seismic data is estimated using a matching point (MP) algorithm having a pre-computed term (e.g., the conjugate of the Gram matrix). A second frequency spectrum of the seismic data is then estimated, using the MP algorithm, an anti-aliasing mask and the estimated first regular frequency spectrum at step 1122. Then, interpolated seismic data is determined at step 1124 by performing an inverse fast Fourier Transform on the second estimated frequency spectrum.

Both of methods 800 and 1100 included discussion of additional post processing steps that should be familiar those of skill in the art. However, as is known to those of skill in the art, there are also pre-processing steps that can be important, and are generally considered to be included in the embodiments discussed herein, and especially with respect to method steps 804, 1104. Interpolation results are known to be sensitive to the quality of input data, so traces used in interpolation should be free from dead or bad traces. The presence of a single dead or bad trace can, under some circumstances, hamper interpolation of all traces generated within that processing window. The distribution of available traces within an interpolation block is also important in the quality of the result. Methods 800 and 1100 according to an embodiment produce good results with input/output traces ratio of about ½ in a 2D case, or about 1/16 in a 5D case, although these embodiments are not limited thereto. According to a further embodiment, distribution of the input traces should be substantially uniform within the interpolation block. Before interpolation according to the embodiment, the input trace data should be verified to ensure a sufficient quantity and that the traces are well distributed. According to a further embodiment, the interpolation windows can be moved, and/or expanded in such regions to minimize edge like artifacts.

Following interpolation, i.e., in method steps 820 and/or 1114, some additional processing can be used to enhance interpolation quality. According to a further embodiment, additional noise removal processes can be implemented, such as prediction error filters or curvelets domain filtering to target noise on interpolated traces before outputting final results.

The methods described herein according to an embodiment for 5D anti-aliasing matching pursuit seismic data interpolation (which collectively refers to both methods 800 and 1100) were tested on a small synthetic dataset containing three linear events, one of which is deliberately steep to test the anti-aliasing capability. FIG. 12A illustrates the original synthetic data, with 32 traces, and FIG. 12B illustrates the same data set with 19 traces randomly deleted. FIG. 12C is an interpolation result using the approach discussed in an article by Naghizadeh, M., et al., 2007, entitled “Multistep Autoregressive Reconstruction (MSAR) of Seismic Records,” and published by Geophysics, 72, V111-E118. The MSAR approach is a combination of an MWNI interpolation method of low frequencies and a prediction error filter applied to high frequencies. The results of interpolation of the data of FIG. 12B using 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment are shown in FIG. 12D.

The interpolation results of both MSAR and 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment are both fairly good, and both methods are able to reconstruct three linear events, including the steeply dipping one. Despite its relative simplicity, however, 5D anti-aliasing matching pursuit seismic data interpolation according to an embodiment compares fairly well against the much more sophisticated MSAR method. The presence of some inconsequential noise can be removed by additional post-interpolation processing as described generally above.

FIGS. 13 and 14 illustrate the results of application of a 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment to field data. Real pre-stack 3D land dataset, sparsely sampled was input to 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment in mid-point and orthogonal offset coordinate space. The bin size was 15 m by 15 m, with in-line offsets ranging between ±1500 m, and cross-line offsets of ±2700 m. The output offset sampling was requested at regular 150 m intervals in the two orthogonal directions, resulting in 21 by 37 traces, i.e. a maximum fold of 777. On average the number of output traces was about 5 times more than the number of input traces. FIG. 13 illustrates a 3D gather before (FIG. 13A) and after (FIG. 13B) application of the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment, wherein traces are ordered by absolute radial offset distance. As can be seen, the input traces were sparsely distributed. As those of skill in the art could appreciate, a conventional 3D interpolation algorithm would not be able to interpolate such an under-sampled dataset. However, the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment is able to fill in all the missing data. Where there are input traces, the interpolated traces converge towards the live traces.

FIGS. 14A and 14B illustrate stacks prior to, and after, use of the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment. The stack results of using the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment shows a significantly better signal-to-noise ratio (SNR) with most of the character being retained from the sparse input data. The noise burst in the shallow data is successfully removed because 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment does not find enough ‘coherency’ in nearby traces to produce a meaningful event.

The 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment described herein is able to match interpolation results from much more complex methods, such as MSAR in synthetic 2D testing. Both methods 800 and 1100 have an inherent simplicity that permits easier implementation. Furthermore, as can be appreciated by those of skill in the art, the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment strikes a balance between the formulation as an optimal weighted normed inverse problem as in MWNI, and the intuitive but non-optimal approach such as ALFT and POCS.

The 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment makes full use of the fact that the Fourier basis functions are the same for all frequency slices. Furthermore, the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment requires no binning, and traces are used at their actual physical coordinates. Still further the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment reproduces original live traces and does not include smoothing, which requires additional pre- or post-processing. This is a crucial advantage for applications that require amplitude information, such as Amplitude Versus Offset (AVO) and Amplitude Versus Azimuth (AVAz) analysis.

FIG. 15 illustrates marine seismic data recording/processing system 1500 suitable for use to implement the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment. Marine seismic data recording/processing system 1500 includes, among other items, server 1501, source/receiver interface 1502, internal data/communications bus (bus) 1504, processor(s) 1508 (those of ordinary skill in the art can appreciate that in modern server systems, parallel processing is becoming increasingly prevalent, and whereas a single processor would have been used in the past to implement many or at least several functions, it is more common currently to have a single dedicated processor for certain functions (e.g., digital signal processors) and therefore could be several processors, acting in serial and/or parallel, as required by the specific application), universal serial bus (USB) port 1510, compact disk (CD)/digital video disk (DVD) read/write (R/W) drive 1512, floppy diskette drive 1514 (though less used currently, many servers still include this device), and data storage unit 1532.

Data storage unit 1532 itself can comprise hard disk drive (HDD) 1516 (these can include conventional magnetic storage media, but, as is becoming increasingly more prevalent, can include flash drive-type mass storage devices 1524, among other types), ROM device(s) 1518 (these can include electrically erasable (EE) programmable ROM (EEPROM) devices, ultra-violet erasable PROM devices (UVPROMs), among other types), and random access memory (RAM) devices 1520. Usable with USB port 1510 is flash drive device 1524, and usable with CD/DVD R/W device 1512 are CD/DVD disks 1534 (which can be both read and write-able). Usable with diskette drive device 1514 are floppy diskettes 1537. Each of the memory storage devices, or the memory storage media (1516, 1518, 1520, 1524, 1534, and 1537, among other types), can contain parts or components, or in its entirety, executable software programming code (software) 1536 that can implement part or all of the portions of the method described herein. Further, processor 1508 itself can contain one or different types of memory storage devices (most probably, but not in a limiting manner, RAM memory storage media 1520) that can store all or some of the components of software 1536.

In addition to the above described components, system 1500 also comprises user console 1535, which can include keyboard 1528, display 1526, and mouse 1530. All of these components are known to those of ordinary skill in the art, and this description includes all known and future variants of these types of devices. Display 1526 can be any type of known display or presentation screen, such as liquid crystal displays (LCDs), light emitting diode displays (LEDs), plasma displays, cathode ray tubes (CRTs), among others. User console 1535 can include one or more user interface mechanisms such as a mouse, keyboard, microphone, touch pad, touch screen, voice-recognition system, among other inter-active inter-communicative devices.

User console 1535, and its components if separately provided, interface with server 1501 via server input/output (I/O) interface 1522, which can be an RS232, Ethernet, USB or other type of communications port, or can include all or some of these, and further includes any other type of communications means, presently known or further developed. System 1500 can further include communications satellite/global positioning system (GPS) transceiver device 1538 (to receive signals from GPS satellites 1548), to which is electrically connected at least one antenna 240 (according to an embodiment, there would be at least one GPS receive-only antenna, and at least one separate satellite bi-directional communications antenna). System 1500 can access internet 1542, either through a hard wired connection, via I/O interface 1522 directly, or wirelessly via antenna 124, and transceiver 1538.

Server 1501 can be coupled to other computing devices, such as those that operate or control the equipment of ship 2, via one or more networks. Server 1501 may be part of a larger network configuration as in a global area network (GAN) (e.g., internet 1542), which ultimately allows connection to various landlines.

According to a further embodiment, marine seismic data recording/processing system 1500, being designed for use in seismic exploration, will interface with one or more sources 118 and one or more receivers 120 These, as previously described, are attached to streamers 116 to which are also attached birds 122 that are useful to maintain positioning. As further previously discussed, sources 118 and receivers 120 can communicate with server 1501 either through an electrical cable that is part of streamer 116, or via a wireless system that can communicate via antenna 124 and transceiver 1538 (collectively described as communications conduit 1546).

According to further embodiments, user console 1535 provides a means for personnel to enter commands and configuration into marine seismic data recording/processing system 126 (e.g., via a keyboard, buttons, switches, touch screen and/or joy stick). Display device 1526 can be used to show: streamer 116 position; visual representations of acquired data; source 118 and receiver 120 status information; survey information; and other information important to the seismic data acquisition process. Source and receiver interface unit 1502 can receive the hydrophone seismic data from receiver 120 though streamer communication conduit 1546 (discussed above) that can be part of streamer 116, as well as streamer 116 position information from birds 122; the link is bi-directional so that commands can also be sent to birds 122 to maintain proper streamer positioning. Source and receiver interface unit 1502 can also communicate bi-directionally with sources 118 through the streamer communication conduit 1546 that can be part of streamer 116. Excitation signals, control signals, output signals and status information related to source 118 can be exchanged by streamer communication conduit 1546 between marine seismic data recording/processing system 1500 and source 118.

Bus 1504 allows a data pathway for items such as: the transfer and storage of data that originate from either the source sensors or streamer receivers; for processor 1508 to access stored data contained in data storage unit memory 1532; for processor 1508 to send information for visual display to display 1526; or for the user to send commands to system operating programs/software 1536 that might reside in either the processor 1508 or the source and receiver interface unit 1502.

Marine seismic data recording/processing system 1500 can be used to implement the 5D anti-aliasing seismic data interpolation method according to an embodiment. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein. According to an embodiment, software 1536 for carrying out the above discussed steps can be stored and distributed on multi-media storage devices such as devices 1516, 1518, 1520, 1524, 1534, and/or 1537 (described above) or other form of media capable of portably storing information (e.g., universal serial bus (USB) flash drive 1524). These storage media may be inserted into, and read by, devices such as the CD-ROM drive 1512, disk drives 1514, 1516, among other types of software storage devices.

It should be noted in the embodiments described herein that these techniques can be applied in either an “offline”, e.g., at a land-based data processing center or an “online” manner, i.e., in near real time while on-board the seismic vessel. For example, data processing including interpolation using the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment can occur as the seismic data is recorded on-board the seismic vessel 102. In this case, it is possible for data generated by the 5D anti-aliasing matching pursuit seismic data interpolation method according to an embodiment to be measured as an indication of the quality of the sampling run.

As also will be appreciated by one skilled in the art, the various functional aspects of the embodiments may be embodied in a wireless communication device, a telecommunication network, as a method or in a computer program product. Accordingly, the embodiments may take the form of an entirely hardware embodiment or an embodiment combining hardware and software aspects. Further, the embodiments may take the form of a computer program product stored on a computer-readable storage medium having computer-readable instructions embodied in the medium. Any suitable computer-readable medium may be utilized, including hard disks, CD-ROMs, digital versatile discs (DVDs), optical storage devices, or magnetic storage devices such a floppy disk or magnetic tape. Other non-limiting examples of computer-readable media include flash-type memories or other known types of memories.

Further, those of ordinary skill in the art in the field of the embodiments can appreciate that such functionality can be designed into various types of circuitry, including, but not limited to field programmable gate array structures (FPGAs), application specific integrated circuitry (ASICs), microprocessor based systems, among other types. A detailed discussion of the various types of physical circuit implementations does not substantively aid in an understanding of the embodiments, and as such has been omitted for the dual purposes of brevity and clarity. However, as well known to those of ordinary skill in the art, the systems and methods discussed herein can be implemented as discussed, and can further include programmable devices.

Such programmable devices and/or other types of circuitry as previously discussed can include a processing unit, a system memory, and a system bus that couples various system components including the system memory to the processing unit. The system bus can be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. Furthermore, various types of computer readable media can be used to store programmable instructions. Computer readable media can be any available media that can be accessed by the processing unit. By way of example, and not limitation, computer readable media can comprise computer storage media and communication media. Computer storage media includes volatile and non-volatile as well as removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CDROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by the processing unit. Communication media can embody computer readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and can include any suitable information delivery media.

The system memory can include computer storage media in the form of volatile and/or non-volatile memory such as read only memory (ROM) and/or random access memory (RAM). A basic input/output system (BIOS), containing the basic routines that help to transfer information between elements connected to and between the processor, such as during start-up, can be stored in memory. The memory can also contain data and/or program modules that are immediately accessible to and/or presently being operated on by the processing unit. By way of non-limiting example, the memory can also include an operating system, application programs, other program modules, and program data.

The processor can also include other removable/non-removable and volatile/non-volatile computer storage media. For example, the processor can access a hard disk drive that reads from or writes to non-removable, non-volatile magnetic media, a magnetic disk drive that reads from or writes to a removable, non-volatile magnetic disk, and/or an optical disk drive that reads from or writes to a removable, non-volatile optical disk, such as a CD-ROM or other optical media. Other removable/non-removable, volatile/non-volatile computer storage media that can be used in the operating environment include, but are not limited to, magnetic tape cassettes, flash memory cards, digital versatile disks, digital video tape, solid state RAM, solid state ROM and the like. A hard disk drive can be connected to the system bus through a non-removable memory interface such as an interface, and a magnetic disk drive or optical disk drive can be connected to the system bus by a removable memory interface, such as an interface.

The embodiments discussed herein can also be embodied as computer-readable codes on a computer-readable medium. The computer-readable medium can include a computer-readable recording medium and a computer-readable transmission medium. The computer-readable recording medium is any data storage device that can store data which can be thereafter read by a computer system. Examples of the computer-readable recording medium include read-only memory (ROM), random-access memory (RAM), CD-ROMs and generally optical data storage devices, magnetic tapes, flash drives, and floppy disks. The computer-readable recording medium can also be distributed over network coupled computer systems so that the computer-readable code is stored and executed in a distributed fashion. The computer-readable transmission medium can transmit carrier waves or signals (e.g., wired or wireless data transmission through the Internet). Also, functional programs, codes, and code segments to, when implemented in suitable electronic hardware, accomplish or support exercising certain elements of the appended claims can be readily construed by programmers skilled in the art to which the embodiments pertains.

The disclosed embodiments provide a source array, computer software, and method for 5D anti-aliasing seismic interpolation by matching pursuit. It should be understood that this description is not intended to limit the embodiments. On the contrary, the embodiments are intended to cover alternatives, modifications, and equivalents, which are included in the spirit and scope of the embodiments as defined by the appended claims. Further, in the detailed description of the embodiments, numerous specific details are set forth to provide a comprehensive understanding of the claimed embodiments. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

Although the features and elements of the embodiments are described in the embodiments in particular combinations, each feature or element can be used alone, without the other features and elements of the embodiments, or in various combinations with or without other features and elements disclosed herein.

This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.

The above-described embodiments are intended to be illustrative in all respects, rather than restrictive, of the embodiments. Thus the embodiments are capable of many variations in detailed implementation that can be derived from the description contained herein by a person skilled in the art. No element, act, or instruction used in the description of the present application should be construed as critical or essential to the embodiments unless explicitly described as such. Also, as used herein, the article “a” is intended to include one or more items.

All United States patents and applications, foreign patents, and publications discussed above are hereby incorporated herein by reference in their entireties. 

We claim:
 1. A method for generating 5D interpolated seismic data using a fast matching pursuit (FMP) algorithm, comprising: obtaining a block of Np irregular trace data, and defining N_(k) to be a desired number of regular trace data; estimating a conjugate of a Gram Matrix G* with size N_(p)×N_(k), using said N_(p) block of irregular traces; transforming said obtained irregular trace data from a t-x domain to an f-x domain; estimating a first set of frequency slices f₁ in an f-k domain of the transformed trace data using a fast MP algorithm that uses G*, wherein f₁ is for low frequency values only; estimating an anti-aliasing mask M based on said first set of frequency slices f₁; estimating a second set of frequency slices f₂ using the fast MP algorithm using G* and the anti-aliasing mask M, wherein f₂ includes all frequencies up to and including a Nyquist frequency of the transformed trace data; and determining N_(k) interpolated trace data by reverse transforming f₂ from the f-k domain to the t-x domain using an inverse 5D fast Fourier Transform.
 2. The method according to claim 1, wherein step (a) of obtaining a block of Np irregular trace data comprises: transmitting seismic waves from one or more sources; and receiving said transmitted seismic waves, and storing the same as irregular trace data.
 3. The method according to claim 2, further comprising: defining said stored irregular trace data into a plurality of interpolation windows; and performing a pre-processing step of moving and/or expanding said plurality of interpolation windows to substantially minimize edge-like artifacts.
 4. The method according to claim 1, wherein said step of estimating a conjugate of a Gram Matrix G* with size N_(p)×N_(k), using said N_(p) block of irregular traces comprises: determining a matrix Φ of size N_(p)×N_(k) by N_(k) column vectors φ_(k), wherein φ_(k) is an exponential wave function k evaluated at x_(l), wherein x_(l) is defined as an instantaneous value of said irregular trace data at a particular point in time; determining a conjugate of the matrix Φ, Φ*; and determining G* according to the following expression— G*=Φ*Φ.
 5. The method according to claim 1, further comprising: increasing a signal-to-noise ratio of said interpolated data by removal of noise from said obtained trace data by determining there to be a lack of coherency between nearby traces containing said noise.
 6. The method according to claim 1, wherein said interpolation occurs without binning of trace data, and further wherein said trace data are used at their actual physical coordinates.
 7. The method according to claim 1, wherein said interpolation reproduces original trace data, and further wherein said interpolation does not substantially smooth said original or interpolated trace data.
 8. The method according to claim 1, wherein said step of estimating a first set of frequency slices f₁ using a fast MP algorithm that uses G*, wherein f₁ is for low frequency values only comprises: initializing a counter n; determining a set of dirty frequency components, F_(D)(n), in the f-k Domain, from the regular trace data, in the t-x domain; defining F_(R)(n) as a first set of frequency slices f₁ of a regular spectrum, and setting F_(R)(n)=0; obtaining a maximum value of the frequency components, F_(D-max)(n), and determining its associated location in the set of dirty frequency components; placing F_(D-max)(n) into a regular spectrum set of F_(R)(n) at a location in the regular spectrum set that is the same as the associated location in the set of dirty frequency components; obtaining a column of G*(n) that corresponds to F_(D-max)(n), and multiplying the G*(n) column by F_(D-max)(n) to obtain G*_(mod)(n); incrementing the counter; subtracting G*_(mod)(n) from F_(D)(n−1) to obtain F_(D)(n); and determining whether there are any new F_(D-max)(n) values, and if yes repeating said steps of determining any new F_(D-max)(n) values until there are no new F_(D-max)(n) values.
 9. The method according to claim 8, wherein F_(D)(n) is calculated according to the following expression: F_(D)(n)=Φ×d, wherein d is a Fourier transform value of a frequency slice of said irregular trace data, and F_(D(n)) is in the f-k domain.
 10. The method according to claim 1, wherein said step of estimating an anti-aliasing mask M based on said first set of frequency slices f₁ comprises: defining an initial anti-aliasing mask M₀ and setting a previous anti-aliasing Mask M_(Prev) equal to said initial-anti aliasing mask M₀; determining an initial residual spectrum RS₀ and setting a previous residual spectrum RS_(Prev) equal to said initial residual spectrum RS₀; multiplying said previous anti-aliasing mask M_(Prev) by said previous residual spectrum RS_(Prev) to determine a residual spectrum-mask product RSMP matrix; determining a maximum value, RSMP_(Max) _(—) _(Val), and corresponding position, RSMP_(Max) _(—) _(Pos), within said RSMP matrix, determining a column in the previously determined conjugate of said Gram Matrix G*, G*_column, that corresponds to RSMP_(Max) _(—) _(Pos), and obtaining wavenumbers from said corresponding G* column; determining a new residual spectrum RS_(New), according to the following expression: RS _(New) =RS−[(α)×(RSMP _(Max) _(—) _(Val))×(G*_column)], where α is a user set variable; and determining that the new anti-aliasing mask M is equal to RS_(New) when [(α)×(RSMP_(Max) _(—) _(val))×(G*_column)] is greater than a threshold T, and when [(α)×(RSMP_(Max) _(—) _(Val))×(G*_column)] is less than or equal to the user set threshold T, setting said anti-aliasing Mask M_(Prev) equal to RS_(Prev), setting RS_(Prev) equal to RS_(New), and determining a new residual spectrum RS_(New) by repeating said steps of determining a residual spectrum-mask product, determining maximum values and positions within said RSMP matrix, and determining a corresponding column within said G*, until said newly determined RS_(New) is greater than the threshold T, or the iteration has reached a maximum number of iterations N.
 11. The method according to claim 10, wherein the number of iterations N are set by a user, and wherein said threshold T is determined by randomly selecting K trace data, wherein K is set by a user, and averaging said randomly selected K trace data to determine said threshold T.
 12. The method according to claim 10, wherein G* is calculated by determining a matrix Φ of size N_(p)×N_(k) by N_(k) column vectors φ_(k), wherein φ_(k) is an exponential wave function k evaluated at x_(l), wherein x_(l) is defined as an instantaneous value of said irregular trace data at a particular point in time; determining a conjugate of the matrix Φ, Φ*; and determining G* according to the following expression— G*=Φ*Φ.
 13. The method according to claim 10, wherein the initial residual spectrum is determined according to the following expression: {tilde over (ƒ)}^(m) =Φ*R ^(m)ƒ; where Φ* is the Conjugate of Φ, and said residual spectrum, RS_(M) is equal to R^(m)ƒ.
 14. The method according to claim 1, wherein use of the conjugate of the Gram Matrix G* reduces computational complexity of the interpolation algorithm from N_(p)×N_(k) to N_(k) operations, wherein N_(p) and N_(k) are a number of input and desired output traces, respectively.
 15. A method for generating interpolated seismic data, comprising: estimating a first frequency spectrum for seismic data using a matching point (MP) algorithm having a pre-computed term; estimating a second frequency spectrum of said seismic data, using said MP algorithm, an anti-aliasing mask and said estimated first frequency spectrum; and determining the interpolated seismic data by performing an inverse fast Fourier Transform on said second estimated frequency spectrum.
 16. The method of claim 15, wherein the step of estimating the first frequency spectrum further comprises: estimating a first regular frequency spectrum for irregular seismic data using a fast matching point (MP) algorithm that includes determination of a conjugate of a Gram Matrix as said pre-computed term.
 17. The method of claim 15, wherein the anti-aliasing mask is estimated by: iteratively multiplying an amplitude of a residual spectrum with a previously determined anti-aliasing mask; determining a position of a maximum element in a result of the multiplication; selecting a coefficient of the residual spectrum at this position, and subtracting a corresponding component of the wavenumber from the residual spectrum to estimate a new anti-aliasing mask.
 18. The method of claim 15, wherein the anti-aliasing mask is estimated using a dip filter estimation method.
 19. A system for generating interpolated seismic data, comprising: at least one processor configured to estimate a first frequency spectrum for seismic data using a matching point (MP) algorithm having a pre-computed term, to estimate a second frequency spectrum of said seismic data, using said MP algorithm, an anti-aliasing mask and said estimated first frequency spectrum; and to determining the interpolated seismic data by performing an inverse fast Fourier Transform on said second estimated frequency spectrum.
 20. The method of claim 15, wherein at least one processor is further configured to estimate the first frequency spectrum by estimating a first regular frequency spectrum for irregular seismic data using a fast matching point (MP) algorithm that includes determination of a conjugate of a Gram Matrix as said pre-computed term. 